%create ref and matrix Up,U0,Vp,V0
clc,clear
close all

%������������������������������for remove fringe��������������������������������������������������������������
YearMonthDay = 20190625;
photoStart=607;
%��������������������������������������������������������������������������������������������������������������
folderName=['G:\Databackup\' num2str(YearMonthDay) 'morning\ikon-22834\'];
Save_Folder_Name = ['D:\oriData\' num2str(YearMonthDay) 'New_Method\'];
mkdir(Save_Folder_Name);
photoNamePre='22834test_';
filename_pre='ascOD22834test_';
refimagesIniDirectory="D:\Rb87ProgramBEC3\removeBackground\backgrounddata\refimagesIni_22834new.mat";
transMatAddress='D:\Rb87ProgramBEC3\removeBackground\RemoveBackgroundNewMethods\transMat.mat';
addpath(genpath('D:\Rb87ProgramBEC3\removeBackground\ExperimentBasicFunction2'));
% %for remove fringe%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% photoStart=350;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
waitSec=4;
drawSwitch=1;
trunc=50000;
% trunc=21000;
%for make background library%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
photoStartRef=101;
photoEnd=350;

refNum=250;
bgmask2=ones(1024,1024);
bgmask3=-ones(1024,1024);
bgmask1=zeros(1024,1024);
ROI=[360 300 290 410];
ROI1=[30 50 970 950];
ROIX_vector=ROI(1):(ROI(1)+ROI(3));
ROIY_vector=ROI(2):(ROI(2)+ROI(4));
ROIX_vector1=ROI1(1):(ROI1(1)+ROI1(3));
ROIY_vector1=ROI1(2):(ROI1(2)+ROI1(4));
bgmask3(ROIY_vector,ROIX_vector)=-2;
bgmask1(ROIY_vector1,ROIX_vector1)=1;
bgmask=bgmask1+bgmask2+bgmask3;
point=[ROI(1) ROI(2);ROI(1)+ROI(3) ROI(2);ROI(1)+ROI(3) ROI(2)+ROI(4);ROI(1) ROI(2)+ROI(4)];
% bgmask1=ones(1024,1024);
% [row,col]=meshgrid(1:1024,1:1024);
% ROI1=[30 50 970 950];
% [bgmask2,point]=ROIfunction(zeros(1024),200,2200,[430 475],[-2 1/2]);
% 
% bgmask1(row<ROI1(1) | col<ROI1(2) | row>ROI1(1)+ROI1(3) | col>ROI1(2)+ROI1(4))=0;
% bgmask=bgmask1.*(~bgmask2);
% figure(11);imagesc(bgmask);
bgmaskL=bgmask(:);
%
if exist(refimagesIniDirectory,'file')
    refimagesIniStruct=load(refimagesIniDirectory);
    refimagesIni=refimagesIniStruct.refimagesIni;
else
    refInd=photoStartRef:min((photoStartRef+refNum-1),photoEnd);
    refimagesIni(:,:,mod(refInd,refNum)+1)=refCreate(refInd,'bg',folderName,photoNamePre,1024);
    save(refimagesIniDirectory','refimagesIni');
end
refSize=size(refimagesIni);
refimages=reshape(refimagesIni,refSize(1)*refSize(2),refSize(3));
refimagesLibrary=refimages;
serial=find(bgmaskL==1);
% refimages=refimages(:);
%%%%%%%%%%%%%%%%incremental SVD
% trunc=0.45;
% trunc=14000;
% trunc=13000;
% trunc=18000;
% refimages=rand(20,30);
% [u1,s1,v1]=svd(refimages);
tic
if ~exist(transMatAddress,'file')
    c=refimages(:,1);
    c=c(serial);
    s=norm(c,'fro');
    U0=c/s;
    V0=1;
    Up=1;
    Vp=1;
    Vpi=1;
    for i=2:refNum
        c=refimages(:,i);
        c=c(serial);
        n=find(c==0);
        if ~isempty(n)
            refimagesLibrary(:,i)=[];
            continue;
        end
        [U0,V0,s,Up,Vp,Vpi]=IncrementalSVD(c,s,U0,Up,V0,Vp,Vpi,trunc);
        U=U0*Up;
        V=V0*Vp;
        figure(10);imagesc(V0);
    end
    toc
    %     save(transMatAddress','U','V','U0','V0','s','Up','Vp','Vpi');
elseif exist(transMatAddress,'file')
    transMat=load(transMatAddress);
    U=transMat.U;
    V=transMat.V;
    U0=transMat.U0;
    V0=transMat.V0;
    s=transMat.s;
    Up=transMat.Up;
    Vp=transMat.Vp;
    Vpi=transMat.Vpi;
    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Remove Background%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
photoInd=photoStart;
while(true)
    savefilename= photoName(folderName,photoNamePre,photoInd,'od.mat');
    if  exist(savefilename,'file')
        photoInd=photoInd+1;disp(['file' photoNamePre num2str(photoInd) 'od.mat' 'exist,continue...']);continue;
    elseif exist(photoName(folderName,photoNamePre,photoInd,'nolight.asc'),'file')
        if exist(photoName(folderName,photoNamePre,photoInd,'nolight.asc'),'file')
            tic
            files_if=[photoInd photoInd];
            [image]=drawpicture(folderName,filename_pre,files_if,'1picture');
            if drawSwitch==1
                figure(1);imagesc(image(1).data,[-0.1,0.5]);title(num2str(photoInd));
                hold on;plot([point(:,1);point(1,1)],[point(:,2);point(1,2)],'r');
            end
            [imSubSize,numZeros]=ifbadpicture(image(1).data,150,150,800,800,150);
            if 3*numZeros>imSubSize(1)*imSubSize(2)
                photoInd=photoInd+1;
                odimage=image.data;
                save(savefilename,'odimage');
                continue;
            end
            [bgRegionStdOri,bgRegionMeanOri]=ODimageMeanStd(image(1).data);
            
            absimages=loadPhoto(photoName(folderName,photoNamePre,photoInd,'.asc'));
            refimagesNew=loadPhoto(photoName(folderName,photoNamePre,photoInd,'bg.asc'));
            refimagesIni(:,:,mod(photoInd,refNum)+1)=refimagesNew;
            nolimage=loadPhoto(photoName(folderName,photoNamePre,photoInd,'nolight.asc'));
            c1=reshape(refimagesNew,refSize(1)*refSize(2),1);
            refimagesLibrary1=[refimagesLibrary c1];
            refimagesLibrary=refimagesLibrary1;
            absIm=reshape(absimages,refSize(1)*refSize(2),1);
            c=c1(serial);
            absIm=absIm(serial);
            
            if InnerProduct(Up(:,1),Up(:,end))>trunc
                [Up]=ModifiedGS(Up);
            end
            [U0,V0,s,Up,Vp,Vpi]=IncrementalSVD(c,s,U0,Up,V0,Vp,Vpi,trunc);
            %             U=U0*Up;
            %             numU=length(find(abs((pinv(U0)-U0'))<1e-15));
            %             disp(['numU=' num2str(numU)]);
            %             V=V0*Vp;
            %             figure(10);imagesc(V0);
            %             refCalCoeff=pinv(V0')*pinv(Vp')*pinv(s)*pinv(Up)*pinv(U0)*absIm;
            refCalCoeff=(V0)*(Vp)*pinv(s)*(Up')*(U0')*absIm;
            refimageCal1=refimagesLibrary*refCalCoeff;
            refimageCal=reshape(refimageCal1,refSize(1),refSize(2));
            refimageCalMNolight=((refimageCal-nolimage)+abs(refimageCal-nolimage))/2+1;
            absimageMNolight=((absimages-nolimage)+abs(absimages-nolimage))/2+1;
            odimage=log(  refimageCalMNolight./absimageMNolight  );
            
            if drawSwitch==1
                figure(2);
                imagesc(odimage,[-0.1,0.5]);
                title(photoName(folderName,photoNamePre,photoInd,'.asc'));
                hold on;plot([point(:,1);point(1,1)],[point(:,2);point(1,2)],'r');%%%%%%%%%%%
                %             rectangle('position',ROI,'Curvature',[1,0],'EdgeColor',[1 0 0.5]);
            end
            [bgRegionStdUpdate,bgRegionMeanUpdate]=ODimageMeanStd(odimage);
            disp(['stdOri=' num2str(bgRegionStdOri) ',meanOri=' num2str(bgRegionMeanOri)]);
            disp(['stdUpdate=' num2str(bgRegionStdUpdate) ',meanUpdate=' num2str(bgRegionMeanUpdate)]);
            
            savefilename= photoName(Save_Folder_Name,photoNamePre,photoInd,'od.mat');
            save(savefilename,'odimage');
            toc
%                         save(refimagesIniDirectory,'refimagesIni');
            
            photoInd=photoInd+1;
        else
            pause(waitSec);
        end
    end
end
    


    
